Self-consistent orbital evolution of a particle around a Schwarzschild black hole 
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The motion of a charged particle is influenced by the self-force arising from the particle's interaction with 
its own field. In a curved spacetime, this self-force depends on the entire past history of the particle and is 
difficult to evaluate. As a result, all existing self-force evaluations in curved spacetime are for particles moving 
along a fixed trajectory. Here, for the first time, we overcome this long-standing limitation and present fully 
self-consistent orbits and waveforms of a scalar charged particle around a Schwarzschild black hole. 
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In spite of the impressive progress made towards tackling 
the two-body problem in general relativity 1 1 1, there remains 
an important regime that appears to be intractable by the meth- 
ods of numerical relativity. When the system consists of a 
massive black hole (M > lO'^Mg) and a stellar mass compan- 
ion (m ^ IQMq), the disparity of length scales characterizing 
this black hole binary proves to be a significant roadblock for 
existing numerical relativity codes. 

The dynamics of such a binary is intuitively simple: A slow 
adiabatic inspiral of the small black hole towards the bigger 
one is followed by an abrupt plunge towards the latter's event 
horizon. However, for the purposes of gravitational wave as- 
tronomy, this qualitative picture is inadequate. The ubiquity of 
supermassive black holes residing in galactic centers has made 
extreme-mass-ratio inspirals (EMRIs) one of the more promi- 
nent predicted sources of low-frequency gravitational waves 
for future space-based missions ||2l. The science we will gain 
from these sources — among them precision tests of general 
relativity in the strong-field regime [3J and a better census of 
black hole populations l|2l |4l — rests on our ability to model 
them to an exquisite degree of accuracy. Specifically, we wish 
to be able to track the phase of their gravitational waveforms 
throughout the long inspiral. 

In the self-force approach to modeling EMRIs, one ignores 
the internal dynamics of the smaller black hole and treats it 
as a massive particle that distorts the spacetime geometry of 
the bigger partner. An EMRI is then equivalent to a charged 
particle moving in a black hole spacetime. But for this ap- 
proach to suffice, the motion of the particle and the resulting 
waveform need to incorporate self-force effects arising from 
the interaction of the particle with its own field. 

Evaluating the self-force is a difficult, though by now well 
understood, process lH] |6l . In a curved spacetime, the field 
generated by a particle at one time backscatters off the cur- 
vature and interacts with the particle at a much later time. 
Consequently, the self-force at any given instant depends on 
the particle's entire past history fT|. This restricts the use- 
fulness of purely analytical self-force calculations mainly to 
astrophysically uninteresting cases [Sl. On the other hand, 
the distributional nature of the point source makes numerical 



evaluation of the self-force technically involved. The retarded 
field diverges at the location of the particle, thus requiring a 
delicate regularization to extract the finite self-force |6|. A 
practical scheme for dealing with this difficulty exists; this 
is the spherical-harmonic-based mode-sum method of Barack 
and Ori [9]. This method has been tremendously successful 
for self-force calculations based on a specified particle orbit 
llT0l - [T3l . However, it has not yet been applied to compute the 
self-force based on an evolved orbit. 

A problem that has resisted solution for a long time is the 
computation of self-consistently self-forced orbits and their 
coiTesponding waveforms. These are orbits that reflect the 
true motion of the particle as it is driven by its actual local 
field. In principle, these self-consistent orbits can be obtained 
only by simultaneously solving the equations governing the 
coupled dynamics of the particle and its field. 

This notion of self-consistency is what we wish to high- 
light, particularly because a recent manuscript |[T4l also re- 
ports on self-forced orbits, though not of the sort we present 
here. In that work, the applied self-force at some instant is not 
what arises from the actual field at that same instant. Instead, 
this applied self-force is what would have resulted if the par- 
ticle were moving for all eternity along the geodesic that only 
instantaneously matches the tiTie orbit; it is the "geodesic" 
self-force and not the self-consistent self-force. It is quite 
likely that the error incurred by this assumption becomes neg- 
ligible in the adiabatic limit, for which the particle stays close 
to the instantaneous geodesic for sufficiently long times. We 
emphasize, however, that this is presently just an expectation 
rather than a demonstrated fact. Rigorously assessing its va- 
lidity requires comparison with fully self-consistent orbits. 

In this Letter, we present for the first time such fully self- 
consistent orbits and waveforms, albeit for a radiating scalar 
charge in the Schwarzschild spacetime. An example is dis- 
played in Figs. [T]and[2] 

Effective source approach:- A novel strategy for self-force 
calculations was proposed 1 15] to address the difficulties aris- 
ing from a 5-function source. Its core idea is to refrain from 
solving the retarded field altogether and to work instead with 
the equation for a field O*^ from which the self-force is readily 
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FIG. 1. (Color online) Orbits for neutral and charged particles start- 
ing at p = 7.2, e = 0.5. The orbital evolution is started close to apas- 
tron (at ? = /i ) and the dots represent events at times ti = 400M, t2 = 
600M,/3 = 1100M,?4 = 1300M,f5 = 1800M and t(, = 2043. 8M, the 
instant of plunge.. The coordinates {x,y} = {r cos (j>,r sin (j)} are 
Cartesian coordinates in the equatorial plane. 



computed: fa = ^V^O^, where q is the charge of the particle. 
The effective source 5'(x"|z"(t), m"(t)) for is C** (contin- 
uous but not differentiable) by construction at the location of 
the particle, in contrast to the traditional 5-function source for 
point particles. Like the 5-function source, S depends on the 
particle's position z"(t) and four-velocity m"(t). (A simi- 
lar approach is independently being pursued in Ref. |fT6ll . the 
main difference there being that it uses a mode decomposition 
in the azimuthal direction). The strategy rests on Detweiler 
and Whiting's insight [17| that there exists a smooth solution 
to the vacuum field equation to which the self-force can be 
fully attributed. This solution is just the difference between 
the retarded field <I>"^' and a locally constructible singular field 
which is the curved spacetime analogue of a "Coulomb" 
field that does not contribute to the self-force. Our approxima- 
tion to the regular field differs from the smooth Detweiler- 
Whiting solution by terms that scale as O(p^) as p — >^ 0, where 
p is some appropriate measure of distance from the particle. 
It is thus only a approximation to the Detweiler- Whiting 
vacuum solution, but it nevertheless gives the same self-force. 
The limited differentiability of <I>^ comes from the inability to 
write down an explicit expression for the full singular field 
from which the effective source is constructed. Generally, 
only approximate expressions for <I>^ are known|ir|. The 
construction of our expression for S is described in detail else- 
where [18!. 

With the effective source at hand, one needs to solve the 
following system of equations: 



dx 

dm 

dx 



:5(x|z(t),m(t)) 



mix) ^ 



(1) 
(2) 

(3) 



where m{x) is the rest mass of the particle. Quinn ||T9l found 
that the rest mass is dynamically modified by the component 
of the self-force tangent to the four- velocity; this is reflected 
in Eq. (|3]l. hi all our simulations, we take the initial rest mass 
m{x ~ 0) to be M. Because of the way S is constructed, <I>^ 
is equal to the retarded field in the wavezone. Thus, by solv- 
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FIG. 2. (Color online) Waveforms from self-consistently evolved 
orbits as detected by an observer located in the orbital plane at 



ing the system of equations above, one obtains not only self- 
forced orbits but their corresponding waveforms as well (see 
Figs. [T]and|2|. 

We recently developed code that solves Eq. ([T]i for a spec- 
ified geodesic in the Schwarzschild spacetime |20'| . Compar- 
ing with Ref. 1 12 1 we find that our main source of error is high 
frequency noise due to nonsmoothness of the effective source 
in the vicinity of the worldline. Most of the time the ampli- 
tude of this noise is small but it reaches a peak of about 2% of 
the value of the self-force at periapsis. 

We then evolve self-consistent orbits by supplementing the 
scalar field evolution with an orbit integrator. Together they 
allow solving Eqs. ([T]l, (|2]l and (j3]l simultaneously. We deal 
with the particle motion in two ways: first by straightforward 
integration of Eq. (|2]l and second by adopting the osculating 
orbits framework described in Ref. f2T|. The first method is 
more general in that it allows us to track the motion of the par- 
ticle all the way to the event horizon. On the other hand, the 
second method, which works only for bound orbits, allows us 
to more readily identify aspects of the evolution that would be 
completely missed by methods relying on flux-averaging and 
balance arguments. For the regimes in which both methods 
are valid, we find the resulting orbits in excellent agreement. 

Self-consistent orbits:- The spherical symmetry of the 
Schwarzschild geometry implies that test-particle orbits may 
always be described by motion in the Q — n/2 plane. The or- 
bits are then characterized by conserved quantities E and L, 
the particle's energy and angular momentum per unit mass, 
respectively. Bound orbits are those for which E < \ and L > 
2y/3M. These possess two radial turning points r± (r_ < r+), 
the periapsis and apoapsis. Following Refs. 112111221 . these or- 
bits can be parametrized in terms of a dimensionless semilatus 
rectum p and eccentricity e, such that r± = pM/{l =F^)- This 
p-e parametrization is geometrically informative: p is a mea- 
sure of the size of the orbit, while e is a measure of deviation 
from circularity. We note, however, that it is meaningful only 
for the space of bound orbits for which {E < 1,L > 2y/3M} 
is mapped onto {0<e< l,p > 6 + 2e}. The separatrix 
p — 6 + 2e corresponds to unstable circular orbits and repre- 
sents the boundary in p-e space separating bound from plung- 
ing orbits. 
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FIG. 3. (Color online) Orbital evolution in p-e space, for an orbit 
that begins with p = 7.2,e = 0.5. Each oscillation in these tracks 
corresponds to one full radial cycle. 




1000 
Time (M) 



2000 



FIG. 4. (Color online) Positional element w interpreted as a measure 
of the periapsis shift relative to its geodesic value. The curves end 
when the particle crosses the separatrix. 



In this parametrization, orbits are described by 

Mp 



d(j) 
dt 



1- 



r{t) 



2Mr' 

r-2M 



1 +ecos{x~w) 



(4) 



[p — 2 — 2ecos(;i; — w)] [1 +ecos{x — w)]^ 



(5) 



where r = r{t), /{t) = ^, p = p{t), e = e{t), % = %(t) and 
w = w(f ) are functions of the KeiT-Schild time coordinate t 
and where x(t) monotonically increases with t. For geodesic 
motion, pif), e(f), and w(f) remain constant; they do, how- 
ever, evolve under the influence of the self-force. 

The orbit we consider here starts at p^ — 7.2, eo — 0.5. In 
Fig.[T] we display an evolved orbit with dimensionless charge 
q := qjM = 1 /32, alongside a test-particle orbit for reference. 
Certain reference points along the orbit are identified to ease 
comparison with our other plots. As initial data, we choose 
<p^(^t ^0) =0 and <i)^(f = 0) = and set the particle initiafly 
moving along the geodesic specified by po = 7.2, eo — 0.5. 
This choice results in a burst of junk radiation, which con- 
taminates the computed self-force at early times but eventu- 
ally goes off the grid. By around t = 200M, the evolved field 
settles down to give the appropriate geodesic self-force, as 
seen in Fig. |5] At f = 400M, when the particle is very close to 
apoapsis, the computed self-force is allowed to act on the par- 
ticle, and the system of equations ([T]l, (|2]), and ([3]) is evolved 
simultaneously for all subsequent times. For this particular 
case (q — 1/32), the particle makes approximately 16 revolu- 
tions (^ 4 full radial cycles) before reaching the horizon. 

Self-forced orbits can also be tracked in p-e space. In 
Fig. [3] we observe that an oscillating and secularly increasing 
eccentricity accompanies the monotonic decrease in p. This 
secular increase is a generic feature of strong-field orbits un- 
der the influence of radiation reaction [22J. The eccentricity 
oscillations, on the other hand, are a new feature of self-forced 
orbits not seen by flux-averaged models. They are due to the 
intrinsic periodicity in the local self-force that goes with the 
(quasi)periodic motion of the particle around the black hole. 
Indeed, it is easy to determine that E = —af^ and L = a^^. A 
decrease in E, while keeping L constant, leads to a decrease 



in e, while a decrease in L, keeping E constant, tends to in- 
crease e. The self-force always decreases both E and \L\, but 
it does so at different rates depending on where the particle 
is. This competition between periodically varying loss rates is 
what leads to the oscillatory behavior in e. 

Since E and L scale as q^, we expect that, starting from the 
same initial conditions, the time it takes a particle to reach the 
separatrix (equivalently, the number of radial cycles) should 
scale approximately as l/q^. This is confirmed by our results. 
For q ^ 1/8, the particle crosses the separatrix and plunges 
before completing one radial cycle. For q = 1/16, 1 /32, and 
1/64, the particle plunges after one, four, and 16 full radial 
cycles, respectively. 

It is difficult to disentangle dissipative and conservative ef- 
fects of the self-force in a self-consistent evolution, since a 
reference geodesic fT3] |23] |24| or an explicit expression for 
the force |21 1 is not available. Nevertheless, it is clear that the 
p-e tracks do not fully characterize the orbital evolution. The 
osculating elements {p,e} are in one-to-one correspondence 
with {E,L}, whose rates of change are determined only by af^ 
and a^^. The r component of the self-acceleration cannot be 
inferred from the p-e tracks alone, and instead its effect man- 
ifests in the secular change of the positional elements 11211 . an 
example of which is w in Fig. [4] The observed changes in w 
are completely missed by flux-averaged approximations. 

The self-acceleration and mass change along the orbit is 
shown in Fig. 1] Only a"^ and are plotted here; the third 
component a' is easily determined from the orthogonality con- 
dition = 0. Inconsistent initial data contaminate the self- 
force early on, but this radiates away before self-consistent 
evolution starts at f = 400M. The self-force depends most sen- 
sitively on the radial position of the particle, with its strength 
increasing the closer the particle is to the black hole. There is 
a small noticeable change in the extrema of the self-force, but 
it is possible that this is mainly due to the small corresponding 
shifts in the extremal radii of the orbit. These properties likely 
describe the geodesic self-force as well; it will be instructive 
to compare self-consistent and geodesic self-forces. 

Earlier we reported an upper bound of 2% on the error in 
our self-force calculation. While this may appear sizable, we 
reiterate that this is merely an upper bound which is reached 
only for very brief portions of the orbit. Our agreement with 
the results of Ref. 1 12 1 is, in fact, significantly better through- 
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FIG. 5. (Color online) Acceleration along a trajectory that starts at 
p = 1 .2,6 = 0.5 for q = 1/32. The blue curve indicates Ma' , the 
red curve M^a'^ , and the gray curve ^ . The bottom plot shows the 
corresponding radial position of the particle. 



first opportunity for these proposals to demonstrate that they 
indeed capture all the relevant features of self-consistent or- 
bits and waveforms. 
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out most of the orbit. By performing a higher resolution 
run, we found that the maximum amplitude of the error was 
halved. The corresponding phase at separatrix crossing, how- 
ever, changed only very slightly (« 0.01%). This suggests 
that the noise in the self-force does not significantly affect the 
phase evolution for the length of the runs we consider here. 

Discussion - We summarize by emphasizing a few points. 
First, our time domain 3D code is versatile: It is not lim- 
ited to low eccentricities, equatorial orbits, or even the 
Schwarzschild spacetime. An ongoing challenge is to devise 
more efficient ways to evaluate the complicated expression 
for the effective source in the KeiT spacetime. This would, 
for example, allow one to check if the recently discovered 
Flanagan-Hinderer resonances flS] persist in a self-consistent 
orbit. Second, our approach makes it possible to assess the 
adiabatic argument on which Ref. lfT4ll is based or, for that 
matter, any other proposal for the computation of self-forced 
orbits. Third, our code readily gives self-forced waveforms at 
(see Fig. |2|. These waveforms did not require any post- 
processing after the computation of the orbit; instead, both the 
orbit and waveform are calculated simultaneously. Finally, it 
is possible to generalize our approach to the more important 
gravitational case. But in that context, we stress that care is 
needed in handling delicate gauge conditions [26 j and possi- 
ble instabilities that may be brought forth by the nonradiative 
low multipoles of the metric perturbation fTT\. This represents 
the next major phase of development for the effective source 
program. 

As is to be expected from a 3D code, the computed self- 
force along an orbit is limited in accuracy compared to other 
methods. Further development will be required to improve 
on this with limited computational resources. Moreover, the 
code is too slow for the task of mass-producing waveforms 
that will sufficiently sample the entire EMRI parameter space. 
However, we emphasize that this is not the objective of our 
approach. The true value of our work lies in its ability to val- 
idate assumptions and predictions arising from all other (pre- 
sumably faster) approximate methods. Our results provide the 
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